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Abstract 

We investigate a Gaussian quadrature rule and the corresponding 
orthogonal polynomials for the oscillatory weight function e luJX on the 
interval [—1,1]. We show that such a rule attains high asymptotic 
order, in the sense that the quadrature error quickly decreases as a 
function of the frequency uj. However, accuracy is maintained for all 
values of uj and in particular the rule elegantly reduces to the classical 
Gauss-Legendre rule as uj — > 0. The construction of such rules is 
briefly discussed, and though not all orthogonal polynomials exist, it 
is demonstrated numerically that rules with an even number of points 
are always well defined. We show that these rules are optimal both in 
terms of asymptotic order as well as in terms of polynomial order. 
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Highly oscillatory integrals are ubiquitous in applied mathematics. The 
numerical evaluation of such integrals has been a rich field of study in the last 
decade, following a series of hallmark papers by Iserles & N0rsett [T2l 
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HD1 [TT] . Prompted by the unexpected succes of modified Magnus expansions 
for oscillatory differential equations [9j, these authors introduced asymptotic 
and Filon-type methods as a generalisation of the classical Filon method 
[3] . Other innovations in this field include the numerical method of steepest 
descent [6], based on numerical evaluation of steepest descent integrals, and 
Levin-type methods [13 [13 [19] , based on solving an associated differential 
equation. A common property of all these methods is that high asymptotic 
order can be attained, by which we mean an error that behaves like O(oj~ a ), 
lo — > oo, for a positive a. Unlike classical asymptotic expansions, whose error 
behaves similarly as lo — > oo, these methods are in principle not limited in 
accuracy for a fixed lo. 

This work concerns an asymptotically optimal oscillatory quadrature rule 
that is valid for all frequencies. The discussion will focus on the Fourier- 
integral 

I[f]:= j\x)e^dx, (1) 

where we shall consider lo ranging from to oo. Of the above mentioned 
quadrature rules, the numerical method of steepest descent can be consid- 
ered to be asymptotically optimal; it delivers an error of size O{L0~ 2n ~ l ), 
when using In quadrature points in the complex plane for this integral. 
Filon-type methods are based on interpolation of the amplitude function /. 
By choosing interpolation points that scale towards the endpoints like 1/lo, 
they deliver O{oo~ n ~ l ) error with the same number of points [lOj. In [5] it 
was pointed out that by choosing the interpolation points for the Filon-type 
method to be precisely the quadrature points in the complex plane of the 
numerical steepest descent method, the so-called superinterpolation points, 
also here an error of 0(uj~ 2n ^ 1 ) can actually be attained with only 2n points. 

Next, one could ask, how well these methods fare when lo is small. Unfor- 
tunately, the superinterpolation points are unbounded in the limit lo — > 0. 
As such, they are not a good choice for small lo. For lo = the Gauss- 
Legendre points are optimal in the sense that they integrate exactly poly- 
nomials of degree 2n — 1 using n points. As for Filon-type methods, the 
interpolation points can be chosen such that they approach the Legendre 
points asw-fO [10]. However, it is not clear what is an optimal way to do 
so. In the context of the so-called exponential fitting methods, such strate- 
gies have been described with relatively few quadrature points [13] , or with 
an heuristic dependence on lo [2]. However, in both cases, using 2n points 
still yields an error of size ©(u; - " -1 ), which is not optimal. 

In a recent work, p], it was shown that Gaussian quadrature rules can 



2 



be constructed for certain oscillatory integral transforms, i.e., integrals over 
unbounded domains. This goes against common wisdom, saying that Gaus- 
sian rules can be found for positive weight functions only. A more correct 
statement is that existence and uniqueness proofs, as well as many construc- 
tion methods, rely on positive weight functions. For integral transforms of 
the form 



it was shown in pQ that applying the Gauss-Laguerre rule along the path 
of steepest descent yields a Gaussian rule, which is exact for / being a 
polynomial of degree < 2n — 1. Moreover, a clear connection was shown 
between polynomial accuracy and asymptotic accuracy for such integrals, 
where higher polynomial accuracy means higher asymptotic accuracy. Thus, 
the Gaussian rule for this integral transform attains an error 0(c<j~ 2 " -1 ). 
Similar results were shown for other oscillators, all resulting in rules with 
nodes in the complex plane. 

The topic of the current work is a further investigation of Gaussian rules 
for oscillatory weight functions, but now on a bounded domain. This, we 
will see, is a more difficult case, since it is no longer possible to remove the 
dependency on u; by a simple rescaling as in the unbounded case. For integral 
([!]), the oscillatory part e luJX is the weight function, and one thus seeks rules 
that integrate polynomials up to degree 2n — 1 exactly. This endeavour 
poses several problems. The rules depend on w in a non-trivial way, and 
must be computed numerically. One cannot guarantee existence of the rules, 
though numerical evidence indicates that all rules with an even number of 
points exist. On the other hand, the resulting rules reduce elegantly to 
Gauss-Legendre by construction in the limit uj — > 0. Moreover, it will be 
proved, under mild assumptions, that the quadrature points tend to the 
superinterpolation points in the high-frequency limit uj — > oo. This implies 
that the method yields optimal asymptotic order in addition to optimal 
polynomial order. These observations are the basis of the statement that 
this Gaussian rule is truly optimal for oscillatory integrals of the form ([I]) 
throughout the frequency regime. 

This paper is built up as follows. In ^TJ we recall some preliminary facts 
regarding Gaussian quadrature and highly oscillatory quadrature methods. 
In analytic expressions for the cases n = 1 and n = 2, as well as several 
numerical experiments, shed light on the questions of existence and other 
properties of Gaussian rules for integrals of the form Q. This section is 
concluded with a set of conjectures, most of which are proved in £j3j on the 
properties of the orthogonal polynomials, and Q on the properties of the 
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quadrature rule. 



1 Preliminaries 



1.1 Gaussian quadrature 

An n-point Gaussian quadrature rule {xj,Wj}™ =1 is a quadrature rule which 
has optimal polynomial accuracy, i.e., it is exact for all polynomials of degree 
< 2n - 1, 

/ f(x)h(x)dx = ^2wjf(Xj), /€7>2n-l. 
Ja 3=1 

It is well known that the nodes of a Gaussian rule are precisely the zeros 
of the n-th orthogonal polynomial with respect to the weight function h 
[2]. Classical theory of Gaussian quadrature ensures that such rules exist 
and that they are uniquely defined whenever h is positive. Moreover, the 
positivity of the weight function also ensures that Xj G [a, b], j = 1, . . . , n, 
and that the weights Wj are all positive. 

The monic orthogonal polynomial p n can be computed in various ways, 
for example by the recurrence 

Pk+i{x) = (x- a k )p k {x) - PkPk-i(x), (2) 

with initial values p_i = and po = 1. Defining the pairing 



(/,<?) 



f(x)g(x)h(x)dx, 



which is an inner product if h is positive, the recurrence coefficients are given 
by 

_ (xPk,Pk) o . (Pk,Pk) 
&k / \ -i Pk / \ • 

{Pk,Pk) {Pk-l,Pk-l) 

Alternatively, writing p n (x) = x n + Ylk=o a k xk , the coefficients a^. 
satisfy the linear system 



(3) 

; O-n-l 



( 

Mi 



Mi 

/'2 



\Mn-l Mn ••• M2n-2/ 

where the moments \x m are defined as 

-6 



ai 



\a n -i/ 



( Mn \ 

Mn+l 
\H2n-l) 



(4) 



Mr, 



x m h(x)dx. 
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1.2 Gaussian quadrature and the method of steepest descent 

In the method of steepest descent, an oscillatory integral of the form ([I]) 
is rewritten along the so-called paths of steepest descent. Assuming / is 
analytic, we may deform the path of integration as follows: 



1 ie 



f(x)e lux dx = / /(-l+ifcrMe-Mt / Hl+itu-^e^dt. 

-1 w J UJ J 

(5) 

Applying the Gauss-Laguerre quadrature on the two resulting integrals yields 
an error that decays like O(oj _2n_1 ), when using n points for each integral 
(so 2n points in total) [6]. The Filon-type method, due to Iserles & N0rsett 
[8], is based on polynomial interpolation of the amplitude function /. In [5j 
it was shown that using the complex points obtained by applying Gaussian 
quadrature on the steepest descent integrals in ^ as interpolation points in 
a Filon-type method, also yields an error decay of 0(w _2n_1 ). Following the 
terminology of |5j, these points are referred to as superinterpolation points: 

Definition 1.1 Let {^jjVj}J=i denote the n point Gauss-Laguerre quadra- 
ture rule. The corresponding 2n superinterpolation points are defined as 

-i + ^luli + ^ 



uj I I UJ 



The corresponding quadrature weights for these points are 



en 



UJ i j=l 



Note that Filon-type methods have the advantage that accuracy can be 
increased by adding interpolation points wherever they are needed, while 
maintaining asymptotic accuracy. This can be used to minimise the effect 
of eventual singularities in the complex plane. This leads, of course, to 
different weights. Also note that the superinterpolation points are not well 
defined in the limit uj — > 0. 



2 Polynomials orthogonal to e iUJX on [—1,1], ana- 
lytic examples and numerical experiments 

Considering polynomials orthogonal w.r.t. the weight e tu)X on [—1,1], the 
non-positive weight function does not enable us to use classical existence 
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and uniqueness theory. However, assuming existence for the time being, at 
least for some n and uj, we denote the n-th monic orthogonal polynomial by 
Pn- The polynomial p% is orthogonal in the sense that, 

J p%(x)x j e iuJx dx = 0, j = 0, 1, . . . , n - 1, u > 0. 

The moments can be computed explicitly, 

fC = J x m e iu > x dx = (-l) m (ia;)- 1 - m (r(l + m, -ioj) - T(l + m, iu)), 

where T(z,a) is the incomplete Gamma functions [17J, or by the recurrence 
w _ e iw - {-l) m e- iul m w u _ 2smu 

Mm — : '■ Mm— 1; Mo — ) 

which is obtained through integration by parts. This recurrence should 
be used with great care, since it's forward stable only when m < uj. For 
m > uj, the backward recursion should be used. Assuming existence of 
these polynomials, they can be computed using the linear system for the 
coefficients Q. For the cases n = 1 and n = 2 the polynomials can be 
computed analytically, and this gives some insight into the nature of the 
higher order rules. 

2.1 The case n = 1 

In the case n = 1 the orthogonal polynomial takes the form pf(x) = x + ao, 
where 

ao = — = I • 6 

(Xq t&nuj uj 

Note that the single quadrature point — oo remains on the imaginary axis. 
The rule is undefined when a; is a multiple of tt, except in the limit to — > 
where ao — > 0. This limit corresponds to the Gauss-Legendre rule. 

2.2 The case n = 2 

In the case n = 2, one obtains for p^ix) = x 2 + a\x + ao the coefficients 
2 + 3cj 2 - 2w 4 + (-2 + uj 2 ) cos(2w) - 4w sin(2a;) 



a 



w 2 (-l + 2a; 2 + cos(2w)) 

2t (-2 + 2w 2 + 2cos(2w) + a; sin(2cj)) 
w(-1 + 2cj 2 + cos(2w)) 
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Figure 1: Two Gaussian points for the oscillatory integral start out at the 

Gauss-Legendre points, ±l/\/3, for uj = and follow these curves in the 

complex plane for increasing oj. 
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Figure 2: Reflecting and zooming in on the left curve in Figure [TJ 
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In the low- frequency limit, for u) — >■ 0, we recover the Gauss-Legendre rule 
again. Observe that 

° 1 = -15 W+ 1575 W 

Since the roots x± satisfy — x + — x_ = a\ and x + x_ = ao, it follows that 
x± = i-^ + 0(o;). They are a perturbation of the Gauss-Legendre nodes 
as expected. 

For w — >• oo, we observe that 



2 cos 2 2a; 2 sin 2a; 



ao = -l + 2 ^ + 0(w" 4 ), 



ai = -— - """7" + — + <9(w~ 4 ), 



2i z sin 2a; 2i sin 2 2ta 
leading to x± = ±1 + r. + O(oj~ 2 ). These turn out to be the two-point 



superinterpolation points as defined in Definition 1.1 — recall that 1 is the 
single root of the Laguerre polynomial of degree 1. 

The roots of the polynomials, which are the two quadrature points for 
the Gaussian rule, are given explicitly by 

x± = (uj (— 1 + 2ui 2 + cos 2w)) 1 i (— 2 + 2a; 2 + 2cos2w + a; sin 2a; ) 

± y/-3 + 6a; 2 - 12a; 4 + 4a; 6 + (4 - 6a; 2 ) cos 2a; - cos 4w + 4a; 3 sin 2a; 

Although this expression is rather complicated, one sees that the points 
exist for all uj, unlike in the case n = 1. Fig. [T] shows the curves that 
the two quadrature points trace out in the complex plane as uj increases. 
The qualitative behaviour seen in this figure also appears for higher n. The 
points leave the real line and drift into the complex plane orthogonal to the 
real line. After an excursion into the complex plane the points attract in 
a rather irregular manner towards the two endpoints of the interval. The 
curves appear to have cusps for certain values of us. Regarding the curves 
as parametric curves, it is clear that a cusp can only happen at a singular 
point, i.e., where vanishes. 

In order to compute the recurrence coefficients, defined in ([3]), that would 
enable us to compute p% , we need the quantity 



(P2,P%)= I ( P Z(x)) 2 e^dx 



16 (2a; 3 cos(o;) + a; 2 (-3 + a; 2 ) sin(a;) + sin(a;) 3 ) 



i v " a; 5 (-l + 2a; 2 +cos(2a;)) 



S 



n 2n 3 7r An 5 n 6n In 

Figure 3: The absolute value of (p%,P2) (continuous) and ^ (dashed) as a 
function of uj. 



For certain values of uj, this expression vanishes and the recurrence coefficient 
«2 does not exist. The first such value occurs near uj = 2tt, and this value 
can be numerically computed to be uj\ = 5.92966 . . .. The values of uj for 
which ,p% ) = are difficult to compute explicitly, but one can give some 
information for large values of uj. Dividing by uj throughout, we obtain 
that (p2 , p% ) = is equivalent to 

2cos(w) / 3 \ . . sin(cj) 3 



Hence, to leading order we have the zeros of sinw as solutions, so uj^ = 
kir + 0(k^ 1 ), k — > oo. The first correction to this estimation gives 

2 

uj k = k-K - — + 0(k~ 3 ) 
kir 

This is illustrated in Figure [3J where also \ is plotted. The figure shows 
that the zeros of {P21P2) appear to coincide with the zeros of |^M, and thus 
these values of uj correspond to the the cusps in Fig. [TJ We may conclude 
from this that there is a breakdown of the recurrence ([2]) for countably many 
values of uj, and that these problematic values correspond to cusps in the 
curves Xj(uj). As a result, p^ is undefined for these special values of u, just 
like pf was undefined for values of uj that are exact multiples of tt. However, 
by avoiding the recurrence and by computing with the moments as in eq. 
Q, for example, one can still compute the orthogonal polynomial of degree 
4 for all values of uj. 
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The quadrature method obtained from the superinterpolation points has 
asymptotic order 3, and this method presumably has the same order. Ap- 
plying it to the integral 

J sin(x)e^'dx, 

the error is shown in Fig. [4| We see that the resulting method indeed 
appears to have asymptotic order 3. Note that a Filon-type method with 
two quadrature points is in general only expected to have asymptotic order 
2. This can be achieved by using the endpoints ±1, or any two points that 
move towards ±1 at a rate of 1/u) |1U| . 

2.3 Numerical experiments for n > 2 

For n > 2, obtaining closed forms expressions for zeros of the orthogonal 
polynomials is highly impractical, if at all possible. Turning to numerics, 
the coefficients of the monic polynomial can be computed from the Hankel 
system Q, and the roots can be found numerically. Note that this system 
is likely to be ill-conditioned. The computations in this paper have been 
carried out in high precision in Maple. 

The roots for n = 3 behave in a similar manner as in the case n = 2. 
In Fig. [5j we see that the cusps in this case seem to coincide with the 
cusps in the case n = 2. However, there is an extra root which is always 
on the imaginary axis, a consequence of the symmetry of the polynomials 
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Figure 5: Three Gaussian points for the oscillatory integral start out at the 
Gauss-Legendre points for u = and follow these curves in the complex 
plane for increasing ui (solid lines). The third point is always purely imagi- 
nary and is not shown. The dashed lines are the trajectories of the roots of 
P2 i as m Fig- [j] 
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Figure 6: Reflecting and zooming in on the left curves in Figure [5] near — 1. 
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with respect to the imaginary axis. In Fig[7J one observes that this root is 
indeed undefined for a set of discrete values lj identified as the cusps in the 
case n = 2. 

One can also compute zeros of polynomials of higher order. For n = 16, 
the zeros behave qualitatively in much the same way as the zeros of ^ . The 
cusps in the curves correspond to the same critical values of lj regardless of 
the root, but these values are different from those in the case n = 2. The next 
experiment concerns the behaviour of the nodes for large uj. Figure [9] shows 
the difference between the 8 roots of pf 6 near —1 and the 8 corresponding 
superinterpolation points. It appears that the difference goes like 0(w~ 2 ), 
similar to what was established for the case n = 2. For large values of u, the 
roots of the orthogonal polynomial tend to the superinterpolation points. 

2.4 Conjectures based on observations 

From the above discussion we can list several conjectures: 

1. Rules with an even number of points exist for all lj. 

2. Rules with an odd number of points are not defined for all lj. 

3. In particular, pfn+i ^ s n °t defined for those critical values of lj that 
correspond to the cusps of the roots of P2 n . 

4. The polynomials are symmetric with respect to the imaginary axis. 
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Figure 9: The absolute difference between the first 8 roots of pf e and the 
corresponding superinterpolation nodes, scaled by a; 2 . 
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5. The roots of p% n approach the 2n superinterpolation points at a rate 
of 0{uj~ 2 ) in the high-frequency limit. 

6. The Gaussian rule based on the zeros of -p^n nas asymptotic order 

In this paper we shall not prove the first conjecture, on the existence of even- 
degree polynomials, but assume it to be true. We shall prove conjectures 3 
(Lemma 3.1) and 4 in $3] and conjectures 5 and 6 in jfij 



3 Properties of the orthogonal polynomials 

In this section we set out to describe a number of interesting properties of the 
orthogonal polynomials that are useful later on to explain the observations 
that were made in the previous section. 



3.1 Symmetry 

Symmetries of weight functions and intervals lead to symmetries of the cor- 
responding polynomials. The complex exponential weight function has sym- 
metries with respect to the imaginary axis, in the sense that w{z) = w{—z). 
Note that the point — z is the reflection of z with respect to the imaginary 
axis. This leads to the following. 

Lemma 3.1 Let w(z) be a weight function such that w(z) = w(—z) and T 
be a contour that is invariant under reflection with respect to the imaginary 
axis, i.e., 

z g r =>- -z sr. 

If a unique monic polynomial p n of degree n exists that satisfies the orthog- 
onality conditions 

j p n (z)z k w(z)dz = 0, k = 0, . . . , n — 1, 

then 

Pn (z) = (-1)>^). (7) 
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Proof 1 Using the symmetry of T and of the weight function respectively, 
we find 

p n (z)z k w(z)dz = J p n (—z)(—z) k w(—z)dz 



(-l) fc / p n {-z)z k w(z)dz 



(-1) / Pn(-z)z k w{z)dz = 0. 



r 



Thus, Pn(—z) satisfies the same orthogonality conditions asp n (z). Since the 
latter is unique, we must have that 



Pn(-Z) = Cp n (z), 

for some constant c. Using that p n is monic yields c = (— l) n . 



The weight function e lU3X satisfies the required symmetry of Lemma 3.1 and 
so does the interval [—1,1]. Therefore the polynomials p^ satisfy the sym- 
metry 0. This implies among other things that they are either purely real 
or purely imaginary along the imaginary axis, depending on the parity of n. 

3.2 Derivatives with respect to u 

We gain useful knowledge about the polynomials by differentiating with re- 
spect to u>. Surprisingly, this differentiation yields the orthogonal polynomial 
of smaller degree. 

Theorem 3.2 The derivative of p^ with respect to uj satisfies 

^(x) = -iPnP^x), n = l,2,..., (8) 
where the proportionality constant, which depends on u, 

o _ (PniPn) (q\ 
\Pn-l'Pn-l) 

is one of the recurrence coefficients ^ of the orthogonal polynomials. 
Proof 2 Define the quantities 

F k (cj) = J x k p%(x)e^ x dx. (10) 
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Note that by orthogonality of the polynomials the first n functions vanish 
identically, 

F k (uj) = 0, fc = 0,...,n-l. 
Differentiating with respect to uj yields 

F' k {uj) = J x k ^(x)e iwx dx + i J x k+1 p" n {x)e iu]X dx =: G k {u) + iF k+1 {ui). 

(11) 

Since F' k {ui) = for k = 0, . . . , n — 1 and Fk+i(to) = for k = 0, . . . , n — 2, 
we must have that 

G k {u>) = f x k ^{x)e iuJx dx = 0, k = 0,...,n — 2. 

Because p^ is monic, is a polynomial of degree at most n — 1. By the 
relations above, it satisfies the exact same orthogonality conditions as p^-\- 
Thus, there must be a constant depending on ui, but not on x, such that 



duj 



(x) = Cn{u))p^_ x {x). 



It remains to determine this constant. Setting k = n — 1 in (11) we find 

= G n -iM + %F n (u) = 0. (12) 
Note from the orthogonality of p^ that 

F n (uj) = ^x n p<Z(x)e i "*dx=(p%,p%). 
Similarly, we find that 



G n -i(oj) = c n (uj) I * x n - l p^x)e^dx = c n ( W )(p-_ 1 ,^_ i : 



Combined with (12) this leads to the result. 



3.3 Three-term recurrence relation 

As it turns out, the coefficients a k and P k of the three-term recurrence 
relation ^ can be computed themselves by a recursion. The following 
result follows from taking the derivative with respect to uj of the recurrence 
relation. Note that a k and j3 k are always functions of uj. In the following 
we frequently omit this dependency and we denote the derivative of a k with 
respect to ui simply by a' k and we do so similarly for f3 k . 
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Theorem 3.3 Let ak and [3 k be the oj- dependent recurrence coefficients for 
Pn, defined by Then we have 

Pk+i = (3k~ ia' k 

and 

otk+i = ak- i -p. ■ 

Pfc+i 

Proof 3 Differentiating the recurrence (|2j) with respect to uj yields 

dPk+i i \ i \®Pk i \ i ui \ of ui t \ o ®Pk-l i \ 

—Q-j-iF) = (a;-a! fc ) — (x) - a k p k {x) - fikP^ix) - p k —^—[x). 

Using Thm. \3.2\ and collecting terms leads to 



(-iPk+i + ~foj)Pk( x ) = d x - a k){-iPk) - /3'k)Pk-i( x ) + if3k-iPkPk-2( x )- 

Matching the leading order coefficients on both sides, using that both p k and 

Pjfe-1 are TOOn * c ; we must have —iftk+i + j^j = — if3 k . The recurrence for /3k 
follows from this. Dividing over yields 

p%(x) = (x-a k - {-if3 k T l fl'k)pt-i{x) ~ Pk-iPt- 2 ( x ), 
and the recurrence for ak follows by comparing to the regular recurrence for 

vl, 

Pk( x ) = (x- a fc _i)p^_ 1 (a;) - /3 fe _ip£_ 2 (x). 

Remark 1 In the theory of orthogonal polynomials, these are sometimes 
known as the deformation equations for the recurrence coefficients. General 
expressions of this kind can be obtained whenever a weight function of ex- 
ponential type is perturbed with a parameter, see for instance \% Prop. 2.1] 
and references therein. 



3.4 Trajectories of the roots 

Finally, we intend to show that the trajectories of the roots in the complex 
plane may have cusps. If so, the norm of the corresponding orthogonal 
polynomial vanishes, and as a result the orthogonal polynomial of one degree 
higher does not exist. 

We start by showing the equivalence between the vanishing derivatives 
of the roots x'j and the vanishing derivatives of the polynomial p^- As in 
the previous section we frequently omit the dependency of the roots on oj in 
our notation. 
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Theorem 3.4 Assume that p^ exists and let Xj(oj), j = 1, . . . ,n, denote its 
zeros. If, for a given uj* , we have 



<V) = 0, j = l, 



,n, 



then 



If Pn is uniquely defined, then the converse is also true. 
Proof 4 Writing p% in terms of its factors, we have 



(13) 



Differentiating with respect to ui yields 
dp%(x) 



j24(u) n (x-xj(u)). 



(14) 



1=1 



From this and condition (13) it follows that = 0, i.e., the partial 

derivative of ' p% vanishes identically. By Theorem 3.2 we find that uj* must 



be a root of f3 n (ui), from which we conclude in turn that (p%,p%) mus t vanish 
at uj = ui* . 



Recall the functions F^{oS) from (10) and their derivatives (11). Letting 



k = n — 1, we find from the above and (11) that 



= F n {uj*) = / pf(x)x n e^ x dx = (pf,pf). 



-i 



To prove the converse, we assume (jp^ ,p% ) = 0, which leads using (11) 
again to 



1 dtfUx) 
_! duj 



x k e iuJ * x dx = 0, k = Q,...,n-l. 



This means that 



satisfies the orthogonality conditions of p^. How- 

ever, it is a polynomial of degree n — 1 and since p^ exists uniquely, the 
lower-degree polynomial can only satisfy the above conditions if it vanishes 
identically. 
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If (Pn >Pn ) = f° r some value of n and uj* , then the monic orthogonal 
polynomial p^+i of degree n + 1 does not exist. In fact, the polynomial p^* 
satisfies all of the orthogonality conditions that p^+i should satisfy, since 

J x k pf (x)e iuJx dx = 0, k = 0, . . . , n. 

This explains why the roots of p^ agree with those of p^ in Figure [5] at 
specific values of uj. The third root, as a function of uj, has poles at these 
values of uj. 

4 Properties of the quadrature rule 

Since the weight function considered here is non-positive, only few of the 
classical results for orthogonal polynomials apply. We do not fully settle the 
questions of existence and uniqueness of the polynomials in this paper. Yet, 
several interesting properties of the quadrature rule can be established and 
we start with the asymptotic order. 

4.1 Asymptotic order 

The following theorem establishes a connection between polynomial accu- 
racy and asymptotic accuracy, on the assumption that the quadrature points 
cluster near the endpoints. Note that a Gaussian quadrature rule with an 
odd number of points in our setting does not qualify here, since due to the 
symmetry at least one root is always on the imaginary axis. In the notation 
of the following theorem, a Gaussian rule with an even number of points 
corresponds to N = 4n and M = 2n, and the theorem thus explains the 
error behaviour we observe in Fig. |4} 

Theorem 4.1 Let {xj, be an 2n-point quadrature rule with polyno- 

mial order N, i.e., 

2n 

and assume that N > In, i.e., the rule is at least interpolatory. 

If the nodes xj can be split in two groups {a;j}™ =1 and {x|}™ =1 , such that 
x 1 - = —1 + 0((jJ _1 ) and x? = 1 + Oiuj^ 1 ), j = 1, . . . , n, then, for an analytic 



0,1..., N - 1. 
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function f the error has the asymptotic decay 

2n 



in „i 
3=1 



00 M - u 



where M = [N/2\ . 

Proof 5 From }16\j we have that an analytic function f can be written as 

f(x) = q 2 M-l(x) + R(x)(x + l) M (x - l) M , 

where q2M-i( x ) * s a polynomial of degree 2M — 1 < N — 1, and R(x) is 
analytic. The error of the method is 

2n 

(x)e wx dx 



zn ~i 
J2wjf(xj) - / /(a 
0=1 



2n „! 

= ^2w j R(x j )(x j + l) M { Xj - 1) M - / R(x)(x + l) M {x - l) M e iulx dx. 
3=1 J ~ 1 
By integration by parts the integral in this expression has asymptotic size 
0{uj~ M ~ l ) . The terms in the sum have a factor which is 0(lo~ ai ), since 
Xj = ±1 + It remains to show that Wj = 0(a; _1 ), j = 1, . . . ,2n. 

Applying the method to the j-th Lagrange polynomial lj(x), we have, 

^ lf\x) . ,=i 



03= I l 1 (x)e^dx = -f^ 
J ~ x k=0 



-iu) k +l 



x=-l 



To conclude that Wj = 0(u 1 ), we need that lj k \±l) = O{oo k ). Now assume 
the node Xj is the first member of the group 1: Xj = x\. Writing lj(x) in 
terms of the two groups of nodes we have 



1 n 2 
x — X; rr x — x~ 



2 ■ 

i=2 X l X i i=l X l X i 

The second of these products is clearly bounded with its derivatives, and we 
can, by Leibniz's formula, concentrate on the first. Here, the denominator 
is of order 0{ijj~ n+1 ). Similarly, the numerator is of order 0(uj~ n+1 ), when 
evaluated in x = — 1. Differentiating will lower the degree of the numerator, 

, n n n 

|: no* n (*-**)■ 

i=2 1=2 i=2,i^l 
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Evaluating in x = — 1 gives that the numerator is of order O(oj~ n+2 ). In 
general the k-th derivative of the numerator is of order 0(uj~ n+k+l ), for 
k > \, so 




Similarly one shows that I- (— 1) = 0{uj k ). The argument can be repeated 
for any Xj, regardless of which group it belongs to. 

4.2 Large oo behaviour of 

Next, a rather remarkable fact will be demonstrated, namely that pointwise 
we have that 



P2n( x ) L n (-iuj(x + l))L n (-iuj(x - 1)), uj ->• oo, (15) 



where L n [x) is the classical Laguerre polynomial of degree n. That is, 
the orthogonal polynomial becomes the product of two rescaled classical 
orthogonal polynomials in the limit uj — > oo. This is, in fact, the polynomial 
that vanishes at the superinterpolation points. To prove this we need the 
following intermediate result: 

Lemma 4.2 Consider a vector of n > 1 points in C, xi, . . . ,x n , and a 
sequence of integers a\ < Q2 < • • . < ot n . We construct the generalised 
Vandermonde matrix G = {x^ J }fj = i, 




G = 




a 2 
1 

»2 

2 




\x: 



n 



X' 



.a 2 

n 




Xi Xj ) , 



l<i<j<n 

where L(x±,X2, ■ ■ . , x n ) is a polynomial 



Proof 6 Consider the determinant as a function of x±, . . . , x. 



det(G) = H(x 1 , ...,x n ). 
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By expanding the determinant along any row, it is apparent that H(x\, . . . , x n ) 
is actually a polynomial. Now 

H(xj,x 2 , ■ ■ ■ ,x n ) = 0, j = 2,3, ...,n, 

since this corresponds to a matrix with duplicate rows. This implies that 
x\ — Xj, j = 2, . . . ,n is a factor of H. Similarly, if one considers H in 
terms of the remaining arguments, one sees that xi — xj, i ^ j, are all 
factors of H . The result follows from this. 



The key to showing something like ( 15 ), is to look at p^ evaluated in the 
superinterpolation points. 

Lemma 4.3 Assume p^ n is bounded in to together with all its derivatives. 
Let Xj, j = 1, . . . , 2n, denote the superinterpolation points. Then 

P2n{x 3 ) = 0(uj- n - 1 ), j = l,...,2n. 

Proof 7 The orthogonality condition for p^n ^ s 
-l 



L 



p 2n (x)x : >e iujx dx = 0, j = 0, . . . , 2n - 1. 

l 



Applying the interpolatory quadrature rule based on interpolation at the su- 
perinterpolation points, along with the assumption on the boundedness of 
p% n , gives Th.3.2] 



2/i 



^WiP2n{xi)4 = 0(uj- 2n - 1 ), i = 0,...,2n-l. 
1=1 

Here the weights are given in terms of the Gauss-Laguerre weights rjj, 

rij 

Wj=w j+n = ^-. (16) 
to 

Denoting y = [wiP2 n (x\),W2P2 n {x2), ■ ■ ■ , W2nP2n( x 2n)] T , this is the Vander- 
monde system 

Vy = 0(u- 2n - 1 ), (17) 
where V = }f7=r By Cramer's rule we have 

! _ Adj(V) 
det(V) ' 
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Entries of Adj(V) are computed in terms of determinants: 



where V 3 ' 1 is V with row j and column i deleted. Using this we can find the 
asymptotics of V^ 1 elementwise. First, using the well known formula for 
the Vandermonde determinant, 

det(V) = ] \ ( r, - Xj) = Cco~< n ^ + 0{oj-< n -^- x ), 
l<i,i<2n 

where C / 0, which follows from the fact that the points belong two groups 
ofn points which are w _1 close to either — 1 or 1. Similarly, if we delete one 
row and one column ofV , one of the groups will be missing a point, and V 3)l 
will be a generalised Vandermonde matrix. The form of the determinant is 
given by Lemma 4-2, and from this it follows that 

det(y J >) = O{oj- {n - 1)2 ), 0<i,j< 2n. 

From this it follows that V~ l = 0(lo ti ~ 1 ) element wise. This further gives, 
from (17), that y = 0(uj~ n ~ 2 ), elementwise, and thus, by (16), the desired 
result. 

Finally we can state the correspondence between roots of p% n and super- 
interpolation points, co — > oo. 

Theorem 4.4 Let Xj, j = 1, . . . , 2n denote the 2n superinterpolation points. 
Assume P2 n exists, and that it is bounded with all its derivatives in uj. For 
each yj such that P2n(yj) = 0, j = 1, • • • , 2n, there is a corresponding index 
I such that 

yj = xi + 0(u)~ 2 ), co —> oo. 
Proof 8 Writing p^ n in terms of its factors, we have 

2n 

P2m( x ) = ]J( X - yj)- 



From Lemma 4-3 we have that 

2/i 



Y[(x k - yj ) = 0(co- n - 1 ), k = l,...,2n. 
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Or, in terms of the Gauss-Laguerre points, 



[T(±l + — - %)(Tl + — - Vj) = o^- 1 ), k 1... . 
oj oj 



From this expression, it is clear that for each j = 1 , . . . , 2n we must have 
yj ~ ±1, which gives O(oj~ n ). The only way C(o; _n_1 ) can be attained is if 



yj = l -^ + 0(uj- 2 ) = Xl + 0(uj- 2 ) 
oj 



for some I. 



Corollary 1 Assume p^ n exists, and that it is bounded with all its deriva- 
tives in oj. The 2n point Gaussian rule, with nodes being the zeros of P2 n , 
is of asymptotic order n + 1. 



Proof 9 From Thm. 4-4 we have that all roots behave like ±1 + O(oj ), 
divided into two equally sized groups. The result thus follows from Thm. 



Note that in the last two results we assumed that the polynomial p^ n is 
bounded in oj, along with all its derivatives with respect to x. This, along 
with the assumption that the polynomials exist for all oj, is not proved in 
this paper. 



5 Conclusions and outlook 

We have presented a Gaussian quadrature rule for integrals on [—1, 1] with 
weight function e tux , oj > 0. The associated family of orthogonal polyno- 
mials Pn is non standard because of the oscillatory nature of the weight 
function. It is shown that for some critical values of oj, the orthogonal 
polynomials of odd degree fail to exist, and that phenomenon is related to 
the fact that the bilinear form defined with respect to e wx is not positive 
definite. However, based on numerical experiments we conjecture that all 
polynomials of even degree exist for any value of oj. We give some results on 
this family of orthogonal polynomials, but their global behaviour in terms 
of oj is bound to be very complicated. We note that they should bridge 
between the Legendre case (when oj = 0) and a product of two rotated and 
scaled Laguerre polynomials in the complex plane, as oj — > oo. 
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We also present results connecting the concept of standard polynomial 
accuracy of Gaussian quadrature rules and that of asymptotic order in terms 
of u), which is of great interest in the context of highly oscillatory quadrature. 
These two ideas are related under the assumption that the zeros of p^[ cluster 
near the endpoints, something that is observed numerically in the paper. 
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